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Abstract 

A novel treatment of non-adiabatic couplings is proposed. The derivation 
starts from the long-known, but not well-known, fact that the wave func- 
tion of the complete system of elctrons and nuclei can be written, without 
approximation, as a Born-Oppenheimer-type product of a nuclear wavefunc- 
tion, X(R), and an electronic one, $^j(r), which depends parametrically on 
the nuclear configuration R. From the variational principle we deduce for- 
mally exact equations for $_R(r) and X(R). The algebraic structure of the 
exact nuclear equation coincides with the corresponding one in the adiabatic 
approximation. The electronic equation, however, contains terms not ap- 
pearing in the adiabatic case, which couple the electronic and the nuclear 
wavefunctions and account for the electron-nuclear correlation beyond the 
Born-Oppenheimer level. It is proposed that these terms can be incorporated 
using an optimized local effective potential. 



I. INTRODUCTION 



The Born Oppenheimer (BO) [1,2] or adiabatic [3] separation of the electronic and nuclear 
motion forms the basis of almost any quantum mechanical study of molecules and solids. 
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The method is based on the large difference between the electronic and nuclear masses, which 
implies that the light electrons are able to adapt almost instantaneously to the configuration 
of the much heavier nuclei. Consider a molecule with N e electrons and N nuclei with masses 
M a and charges Z a and denote by r the set of electronic coordinates r — r± . . .r Ne and by 
R the set of nuclear coordinates R = R 1; . . . ~R N . 

The separation is effected by decomposing the complete molecular wavefunction, 
^moi(r,R), in terms of a product of an electronic wavefunction depending parametrically 
on the nuclear positions and a nuclear wavefunction: 

* mol (r,R)^® B R °(r)X BO (R). (1) 

Since the electrons, more or less, cannot see the nuclear motion, the electronic wavefunc- 
tion is chosen to be an eigenfunction of an electronic Hamiltonian derived from the molecular 
one by ignoring the nuclear kinetic energy or, equivalently, the nuclear mass is set to in- 
finity. Finally, following Longuet-Higgins [3] and Messiah [4], the equation for the nuclear 
wavefunction can be obtained variationally. Keeping the electronic wavefunction fixed, one 
requires that the expectation value of the molecular Hamiltonian in terms of the adiabatic 
product wavefuction remains stationary when the nuclear wavefunction is varied. In the re- 
sulting effective nuclear Hamiltonian, in the absence of electronic degeneracies, the nuclei lie 
in a potential which, to high accuracy, coincides with the electronic energy eigenvalue. Each 
electronic level allows for a set of nuclear vibrational and rotational excitations but does not 
couple directly to any nuclear state of this set. For many physical systems and phenomena 
this separation is extremely successful, at least whenever degeneracy or near degeneracy 
of the electronic Hamiltonian is not involved. There are however other phenomena, which 
crucially depend on the coupling between the electronic and the nuclear wavef unctions. For 
these phenomena the adiabatic or BO approximation breaks down. 

One should be careful about the variational derivation of the nuclear BO equation, be- 
cause the ground state of a free molecule is a continuum state (since the complete molecular 
Hamiltonian is translationally invariant). Rigorously then, one should first separate off the 
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centre of mass of the molecule and then perform the adiabatic decomposition of the resulting 
bound wavefuction [5]. This leads to complicated terms in the resulting Hamiltonian. In 
this paper we do not separate off the centre of mass first, because our main aim is to show 
how one may improve upon the BO scheme, keeping at the same time the formulae simple 
and transparent, rather than maintaining complete rigour. One may justify the approach 
by imposing a confining potential that bounds the molecule within a very large but finite 
volume. 

In probability theory, the joint probability distribution p(r, R) of two sets of variables r, R 
can be factorized as the product p(r, R) = p(r\R) p(R) of a conditional probability distribu- 
tion p{r\R) for the (electronic) variables r given R, times a marginal probability distribution 
p(R) for the (nuclear) variables R. Hunter, proposed accordingly that the molecular wave- 
function ^ mo \(r, R), i.e., the joint probability density amplitude of electronic and nuclear 
degrees of freedom could also be decomposed exactly as a product of a conditional probabil- 
ity amplitude &n(r) of electronic variables, given fixed values for the nuclear variables times 
a marginal probability amplitude X(R) for the nuclear degrees of freedom [6]: 

* mA (r,R) = Q R (r)X(R). (2) 

The BO electronic and nuclear wavefunctions ^^(r), X BO (R) provide excellent approxi- 
mations for &n(r) and X(R) but Hunter suggested that for each \l/ mol (r, i?) one may find 
a conditional and a marginal probability density amplitudes, such that the factorization is 
exact. Czub and Wolniewicz discovered that for diatomic molecules, X(R) is a nodeless 
wavefunction and hence the nuclear non-adiabatic potential must necessarily demonstrate 
spikes at the points were the adiabatic wavefunction has nodes [7] . Hunter observed that this 
result is more general [8,9]. The proof is straightforward: Take a molecular wavefunction 
that is to a good approximation given by a single BO product, ty mo \(r, R) = $^(r) X BO (R), 
with the electronic part in the ground state and the nuclear wavefunction in an excited state 
with a node at Rq. Since the electronic adiabatic states form a complete set, the molecular 
wavefunction can be expanded exactly in the series. 
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* mol (r,R) = J2*n°R(r)Xn(R). (3) 

n 

Consequently, the diagonal of the nuclear iV-body density matrix is equal to 

(* mo l(i?)|* mo l(i?)> = E \Xn{R)\ 2 . (4) 

n 

Here and in the following, integration over all electronic degrees of freedom is denoted by the 
bra-ket symbols, while integration over the nuclear degrees of freedom is written as / d 3N R. 
In order for {9moi(R)\^moi(R)) to have a node at Ro, all vibrational states X n (R) must have a 
node at the same R , which cannot happen. As an interesting aside, we note that Hunter's 
argument fails for nodes dictated by the symmetry. For example for identical fermionic 
nuclei, the molecular wavefunction must be antisymmetric with respect to exchange of two 
nuclei, say at positions Ri, R2 and then the diagonal iV-body nuclear density matrix and 
consequently the marginal probability amplitude must have nodes at (Ri,R 2 ) = (R, R). 

So far, the impact of the work on the exact factorization of the molecular wavefunction 
has been limited and these ideas have not led to new practical ways to improve upon the 
adiabatic approximation. Perhaps in the early papers, attention focused too heavily on the 
nodeless character of the marginal probability density amplitude, which admittedly is a non- 
adiabatic coupling but probably not a very important one, since when the BO approximation 
is accurate, the correction does not have a measurably significant effect and when the BO 
approximation breaks down, it is not due to the nodeless property of the wavefunction. More 
importantly and rather curiously, the non-adiabatic equation determining the electronic 
wavefunction $i?(r) was never given. Knowledge of that equation would inform us of the 
non-adiabatic terms that correct the BO electronic wavefunction and hopefully accurate 
approximate treatments of these terms would emerge. This is the goal of our paper. 

There is a representability issue which perhaps was not addressed adequately, namely, 
it is not clear whether the marginal probability distribution always corresponds to a 
meaningfull nuclear wavefunction. For example the decomposition that was proposed 
to demonstrate the factorization (2) with $i?(r) = ^ mo i(r, R)/\J mo\(R)\^ m.o\(R)) and 
X(R) = y/(* mol {R)\* mol {R)) 
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implies that the nuclear wavefunction is always bosonic even when the nuclei are 
fermionic. This problem can be overcome by introducing an /^-dependent phase and writing 

<^mol(i?)> yJ(* m0 l(R)\* m0 l(R)) 



where g(r) is a suitably chosen electronic wavefunction that does not depend on R. In 
this way, the electronic wavefunction $i?(r) (conditional probability amplitude) depends 
on R and for identical nuclei is symmetric with respect to nuclear exchange. The nuclear 
wavefunction X(R) (marginal probability amplitude) has the same symmetry properties as 
$ mo i(r, R) with respect to nuclear exchange. 

Adopting the small potential mentioned earlier to confine the molecule in a very large 
volume, the molecular, electronic and nuclear wavefunctions are all normalized to unity for 
convenience: 

J d m R(V mol (R)\V mol (R)) = 1, (7) 
($ R \$ R ) = 1, (8) 

J d 3N R\X(R)\ 2 = 1. (9) 

The electronic and nuclear wavefunctions have the freedom of an R— dependent phase 
(which, for identical nuclei must be symmetric with respect to exchange) 

&r(t) -> e ld{R) $> R (r) , X(R) -> e~ im X{R) (10) 

since this transformation leaves *& mo i(r, R) invariant. Apart from this freedom of phase the 
components X and are unique. 

II. NON-ADIABATIC ELECTRONIC AND NUCLEAR EQUATIONS 



As the parametric product ansatz (2) does not constitute an approximation, it should be 
possible to improve significantly upon the BO level of accuracy, still using a product wave- 
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function (2), if we optimize the electronic wavefunction variationally, under the normalisation 
constraint (8), rather than choose it a priori to be an eigenstate of the BO electronic Hamil- 
tonian. Naturally, we expect that in the non-adiabatic scheme, the electronic wavefunction 
will be coupled to the nuclear state. 

The molecular Hamiltonian has the form: 

N h 2 X7 2 

H mol = H%° - V ^ (11) 
£i 2M a 1 ; 

where V„ is the Laplacian with respect to R a and H®° is the electronic Hamiltonian: 



H 



= -^Ev^2S--EE F ^n + 2E^f- (12) 

Z '" L j=l z i,j=i ' ij j=la=l\ L 3 "-al z a,0=i n-aP 

The nuclear Coulomb energy (last term) is independent of r and represents a constant shift 
of the electronic energy. 

In order to obtain the equations determining the electronic and the nuclear wavefunc- 
tions, $i?(r) and X(R), we have to vary the expectation value of the Hamiltonian H mo \ 
with respect to $i?(r) and X(R), under the normalisation constraints (8, 9). The latter are 
incorporated through the use of Lagrangian multipliers A'(R), E. Note we have infinitely 
many constraints for the normalisation of the electronic wavefunction, one for each R, and 
we need as many Lagrangian multipliers A'(R): 

-Jd 3N RA'(R) 1) , (13) 

—E J d 3N R(\X(R)\ 2 -1). (14) 

A few words about notation: If the gradient operator V a , or the Laplacian V^, lies within 
parentheses or between bra-kets, as in (<&r| V a $ij), it is implied that it acts only inside the 
bra-ket or parentheses. Note that when we use square brackets, as in [— i%V a ± A a ], the 
action of V a is not confined within the brackets. Direct optimization of the electronic and 
nuclear wavefunctions under the constraints (13, 14) yields the non-adiabatic equations: 



h bo 1 ^ X*(R)(-ihV a X(R)) 



$ R (r) = A(R) 0> R (r) (15) 
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[ - ihVg + A a (R)Y , bo,^ u f Wr|VA! 



.a=l 



a=l 



2M n 



X(i?) =EX(R). 
(16) 



where A(i2) = M(R)/\X(R)\ 2 and A a (R) is the vector-potential-like real quantity 



A Q (i2) = ($ fl |-iWA>. 



To make further progress, we define an electronic energy as 



(17) 



•r 



Q = l 



N HhV a -A a (R)} 2 



2M n 



R 



(18) 



which, by Eq. (17), can be written as 



^(i?)-($ R |^ !<*>*> + E ^- E^~- 



(19) 



After some straight-forward algebra, Eqs. (15) and (16) can be expressed in terms of this 
electronic energy as 



K° + E HWa " Aa{R)? 

a=l 

N 1 



2M n 



' 5 * . 



X*(R) (-ihV a X(R)) 

\*(R)\ 2 



+ A a (R) 



ihV a -A a (R) 



$ R (r) = £(R)$ R (r) (20) 



V HW a + A a (i?)] 2 
^ 2M + ( } 

La=l Z1Vln 



X(R) = EX{R). 



(21) 



While, from the numerical point of view, this form of the equations is clearly not ad- 
vantageous, the new form reveals some other interesting aspects: First of all, Eqs (20) and 
(21) are manifestly form-invariant under the gauge transformation (10). In particular, the 
electronic energy (18) which appears as "eigenvalue" in (20), is gauge invariant. The third 
term on the left-hand side of Eq. (20), while essential for the coupling between X and 
$r, does not contribute to the energy £(R), i.e., when Eq. (21) is multiplied from the left 
with $^(r)* and integrated over r, this term drops out. One may confirm directly that the 
solution of the non-adiabatic equations (15, 16) is exact by showing that if $«(r) and X(R) 
satisfy (15, 16) then their product will be an eigenstate of H mo \. 



H mol $ R (r) X(R) = E $ fl (r) X(R). (22) 

If the phase of the electronic wavefunction $i?(r) depends on r only, and if the dependence 
on R is continuous and differentiable, then the vector potential A a (R) vanishes and (16) 
simplifies to 



£^ + M o^ ) + ^<VA|VA> 



X{R) = EX{R). (23) 



L a=! 2M a * « 1 2M a 

The algebraic structure of the non-adiabatic nuclear equations (16, 21, or 23), coincides 
with the nuclear equation one encounters in the adiabatic approximation. However, in the 
latter, the electronic wavefunction is the BO one, satisfying 

H B ^ BO (r)=£ BO (R)^ BO (r), (24) 

with £ BO (R) = (^°\H^°\^°), while in the non-adiabatic nuclear equations (16, 21, or 
23), the electronic wavefunction $i?(r) is the non-adiabatic one satisfying (15) or (20). 

The variational derivation of Eq. (16) guarantees that even in the approximate adiabatic 
scheme, the lowest energy eigenvalue, denoted by E BO , is rigorously an upper bound to the 
exact ground state energy of the system: 

E < E BO (25) 

since E is the minimum of the molecular Hamiltonian in a Hilbert space which is wider than 
the Hilbert space of the adiabatic states. The inequality holds even when the electronic 
non-adiabatic equation is derived in a constrained way, as in the following section on the 
Optimized Effective Potential Method treatment of the non-adiabatic terms. It is interesting 
to note that 

($ R \H%>\$ R )>£ BO (R) (26) 

since S BO (R) is, for each R, the minimum of the expectation value of H B °. 
A necessary consequence of (25) and (26) is that 
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£ 2M I d3NR \ X ( R )\ 2 ((Va^|VA> - <V a $!°|V a $!°>) < (27) 

where X(R) is the lowest non-adiabatic nuclear state. Observe that all terms in the inequal- 
ity are nonnegative. The inequality indicates that "on average" for the non-adiabatic elec- 
tronic state that corresponds to the lowest nuclear state (V a &n\ V a ®R) < (Vq,^ ! V a $^°). 
So, $r can be expected to change less rapidly than <3>^° as a function of R. It is worthwhile 
to remember that when a nucleus samples a region of a degeneracy or an avoided crossing 
of electronic levels, the concept of the adiabatic surface looses its meaning and an indication 
is that the potential energy depends strongly on the nuclear position. Small changes in R 
result in large changes in the wavefunction and potential. A lot of effort has been invested 
in constructing a different basis in the space of (near) degenerate electronic wavefunctions 
with energies that would behave smoothly close to the (avoided) level crossing. A conve- 
nient choice is the diabatic basis [12], defined as the basis that diagonalises the off-diagonal 
matrix elements ($_R, m | V a $_R in ) [12,13]. It is known however, that rigorously, such a basis 
does not exist in general [13,14]. The non-adiabatic electronic wavefunctions and energies 
may offer a physical solution to this problem. There are cases when the the vector-like- 
potentials A a cannot all be made to vanish (for example when adiabatic electronic levels 
become degenerate at a point R), leading to a pseudo-magnetic field in (16) and a geometric 
phase [10]. An interesting question is whether the geometric phase is a consequence of the 
approximate nature of the adiabatic approximation and consequently whether it would dis- 
appear in the exact non-adiabatic scheme. We believe the geometric phase remains in the 
non-adiabatic scheme, since the non-adiabatic corrections depend on the nuclear mass and 
vanish for M — > oo, whereas the geometric phase is independent of the size of the nuclear 
mass. The geometric phase is then not an artefact of the adiabatic approximation but a 
result of the parametric decomposition (2) of the molecular wavefunction which, by itself, 
is not an approximation. 
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III. APPROXIMATIONS FOR THE ELECTRONIC NON-ADIABATIC STATES 



The exact non-adiabatic electronic equation is naturally very hard to solve and one has 
to resort to approximation techniques. A perturbative treatment of the non-adiabatic terms 
in the electronic equation is straightforward to apply if we keep the nuclear wavefunction 
X(R) fixed. However, when the non-adiabatic terms are large (which is the case of interest 
here), the expansion may converge very slowly, or not converge at all. Another method 
to deal with the non-adiabatic electronic equation, that seems particularly promising, is to 
seek for the local independent potential Vr(t) that best approximates the effect of the non- 
adiabatic terms in (15), in a manner analogous to the Optimised Effective Potential (OEP) 
method [15-17]. 

Remember that the non-adiabatic equation was derived under the constraint of normal- 
isation of the electronic wavefunction (8). This was the most general way of requiring that 
the electronic wavefunction depends parametrically on the nuclear positions. A more re- 
strictive but perhaps more physical way is to impose that the electronic wavefunction must 
be an eigenstate of a Hamiltonian which is the sum of H® and an electronic potential Vr 
that depends on the nuclear positions. The variationally best potential is determined by 
making the expectation value of the molecular Hamiltonian (11) stationary. The optimal 
Vr will in general be an electronic A e -body potential. A further approximation, to make 
the problem tractable, is to restrict to (local) one-body potentials. Hence we optimize the 
electronic wavefunction $r(v) under the condition that it satisfies the equation 

& R (r) = E(R)$ R (r). (28) 

In this way, the expectation value of the molecular Hamiltonian (11) becomes a functional 
of the effective non-adiabatic potential, (H mo \) = E[V R }. 

Varying the effective potential, we obtain the functional derivative of £"[Vr] with respect 
to V R : 
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H 



BO 
R 



MM, 



+ | i i™ 1 -<W)G E |VA) + c.c. (29) 

where, 

G * - \ € U (R)-€(RY (30) 

Summation is over all eigenfunctions of, + Vr, except <3>#. For compactness, we use 
in (29,30) second-quantization notation, p(r) is the electron charge density operator and 
V R = f<FxV R (x)p(x). Obviously, G R \$ R ) = 0. 

The optimized effective potential we seek corresponds to a stationary point of ^[Vr], i.e. 
the functional derivative vanishes 

^1 - (31) 

Solving the OEP equation (31) to obtain the potential can be very expensive. The 
approximation proposed by Krieger, Li and Iafrate [18,19] (KLI) can be used to simplify the 
equation considerably. 

When the electronic state is a Slater determinant of Kohn-Sham spin-orbitals (f Rta , one 
obtains the spin-dependent KLI effective potential: 

2^>K(r) = 2 / rfrV fi , ff (r')ElO)lVyr')| 2 

J 3 

h 2 * h 2 

- E ^TT E V^/MVjy^r) |</4 !CT (r)| V RiCT | Vjy^) + c.c. 

- E W ^f^ ■ ? ^,/(r)V a ^ CT (r) - V*,.|V^>] + c.c. (32) 

This equation can be solved within the ordinary Kohn-Sham self-consistency loop: The 
potential Vr^{y) to be used in the next iteration is obtained by plugging the Kohn-Sham 
orbitals and the Vr >(7 (y) of the previous iteration in the right-hand side of Eq. (32). 

We have solved the non-adiabatic OEP equation for the simplest system that contains 
non-adiabatic couplings, namely the hydrogen ion H^- Using the BO scheme one obtains 
an excellent approximation for the ground state wavefunction, \&(r, R) = $^°(r) X BO (R). 
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To test our method, we used the inverse BO separation, i.e., we wrote, \l/(r, R) = 
0i. BO (R) X lBO ( r ) ? where (Jy^CR) is the conditional probability density amplitude for the 
relative nuclear coordinate, when the electron is frozen at r, and x lBO ( r ) the marginal prob- 
ability density amplitude for the electron, moving in the potential surface of the nuclei. The 
inverse BO separation ignores terms of the order M/m and is hence very poor. Including the 
non-adiabatic correction terms (in this case the non-adiabatic terms can be transformed to 
have the form of a potential and hence the OEP method is exact), we found that the solution 
converged numerically to the correct wavefunction one obtains doing the BO seperation in 
the orthodox way. 

IV. SUMMARY 

We have revisited the adiabatic/BO separation of the electronic and nuclear degrees of 
freedom in molecules, aiming to improve accuracy beyond the adiabatic level, but keeping 
at the same time the language and concepts of the adiabatic approximation which form the 
basis of almost all studies of molecular systems. We have derived a non-adiabatic correction 
for the electronic equation, which couples in a self-consistent way the electronic and nuclear 
wavef unctions. This correction has the form of an effective potential. The implementation 
of this correction to electronic structure codes that already use the OEP method to account 
for electronic exchange and correlation seems particularly promising. 
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